home *** CD-ROM | disk | FTP | other *** search
/ SGI Developer Toolbox 6.1 / SGI Developer Toolbox 6.1 - Disc 4.iso / src / haeberli / mtex / resample.c < prev    next >
C/C++ Source or Header  |  1994-08-01  |  15KB  |  400 lines
  1. /*
  2.  * Copyright 1991, 1992, 1993, 1994, Silicon Graphics, Inc.
  3.  * All Rights Reserved.
  4.  *
  5.  * This is UNPUBLISHED PROPRIETARY SOURCE CODE of Silicon Graphics, Inc.;
  6.  * the contents of this file may not be disclosed to third parties, copied or
  7.  * duplicated in any form, in whole or in part, without the prior written
  8.  * permission of Silicon Graphics, Inc.
  9.  *
  10.  * RESTRICTED RIGHTS LEGEND:
  11.  * Use, duplication or disclosure by the Government is subject to restrictions
  12.  * as set forth in subdivision (c)(1)(ii) of the Rights in Technical Data
  13.  * and Computer Software clause at DFARS 252.227-7013, and/or in similar or
  14.  * successor clauses in the FAR, DOD or NASA FAR Supplement. Unpublished -
  15.  * rights reserved under the Copyright Laws of the United States.
  16.  */
  17. #include <gl/gl.h>
  18. #include <time.h>
  19. #include <math.h>
  20. #include "types.h"
  21. #include "image.h"
  22. #include "colors.h"
  23. #include "other.h"
  24.  
  25. #define EPS .001
  26. Boolean d;
  27. typedef unsigned int unsi;
  28.  
  29. extern IMAGE *iopen(char *,char *,unsi type,unsi dim, unsi xsize, unsi ysize, unsi zsize);
  30. IMAGE *oimage;
  31. kernel_list_t kernel;  /* Global...created by generate_kernels, used by sample_gaussian */
  32.  
  33. /***************************************************************************/
  34.  
  35. void generate_kernels(float st_dev, float threshold)
  36. {
  37. long i,k,region,sgrid,rsize,osize,ksize,ssize,kx,ky,x,y,coeff_count;
  38. float sx,sy,dx,dy,dist,dist_scaled,dev_sq,gauss,sum,inside_sum,scaled_threshold;
  39. float peak_gauss;
  40. float *kernel_tmp,*kc;
  41. long *ks;
  42. short *ko;
  43. Boolean exclusion_ok;
  44.   
  45. /* Make an array of convolution kernel coefficients.  Array will have many    */
  46. /* kernels, one for each of a subpixel grid.  Both the number of coefficients */
  47. /* and the size of the subpixel grid are dynamic based on the standard        */
  48. /* deviation, which is presumably the zoom factor.  Each kernel will have a   */
  49. /* variable # of coefficients...whatever nearby pixels have coefficients over */
  50. /* threshold will be included.                                                */
  51.    
  52.    if (st_dev < .4) st_dev = .4;   /* Don't undersample if magnifying */
  53.    region = (long)(10.0 * st_dev);
  54.    if (region < 7) region = 7;
  55.    sgrid = (long)(16.0 / st_dev);
  56.    if (sgrid < 1) sgrid = 1;
  57.    kernel.subpixel_grid_size = sgrid;
  58.    
  59. /* Dynamically allocate storage for arrays pointed to by kernel structure */
  60.    
  61.    ssize = (SQR(sgrid) + 1) * sizeof(long);   
  62.    rsize = SQR(region) * sizeof(long);
  63.    ksize = SQR(sgrid) * SQR(region) * sizeof(float);
  64.    osize = SQR(sgrid) * SQR(region) * 2 * sizeof(short);
  65.    kernel.starts =   (long *)malloc(ssize);   /* One "start pointer" per subpixel gridpoint, plus extra at end */
  66.    kernel_tmp =     (float *)malloc(rsize);   /* One kernel_tmp entry per nearby pixel                         */
  67.    kernel.offsets = (short *)malloc(osize);   /* Two offsets per coefficient                                   */
  68.    kernel.coeffs =  (float *)malloc(ksize);   /* One coeff per nearby pixel per subpixel gridpoint             */
  69.    ks = kernel.starts;
  70.    kc = kernel.coeffs;
  71.    ko = kernel.offsets;
  72.    if (kernel.coeffs == NULL)
  73.       {
  74.       printf("Unable to allocate %d bytes for resampling kernels\n",ksize);
  75.       exit(1);
  76.       }
  77.    coeff_count = 0;
  78.    for (ky=0; ky<sgrid; ky++)               /* Subpixel offsets...each will be kernel */
  79.    for (kx=0; kx<sgrid; kx++)
  80.       {
  81.       sum = 0.0;
  82.       peak_gauss = 0.0;
  83.       k = 0;
  84.       sx = (float)(region/2) + (float)kx / (float)sgrid;   /* Subpixel center */
  85.       sy = (float)(region/2) + (float)ky / (float)sgrid;
  86.       
  87. /* Calculate gaussian coefficient of distance from x,y to sx,sy */     
  88.       
  89.       for (y=0; y<region; y++)              /* Consider nearby pixel centers */
  90.       for (x=0; x<region; x++)
  91.          {
  92.      dx = sx - (float)x;                /* X distance from pixel center to subpixel center */
  93.      dy = sy - (float)y;                /* Y distance from pixel center to subpixel center */
  94.       dist = sqrtf(SQR(dx) + SQR(dy));   /* Euclidean distance */
  95.      dist_scaled = dist / st_dev;
  96.      dev_sq = SQR(dist_scaled);
  97.          gauss = exp(-0.5 * dev_sq);        /* Gaussian coefficient for this dist */
  98.      kernel_tmp[k++] = gauss;
  99.      sum += gauss;
  100.      peak_gauss = MAX2(gauss,peak_gauss);
  101.      }
  102.      
  103. /* Sum the coefficients which exceed threshold so we can normalize them */    
  104.      
  105.       inside_sum = 0.0;
  106.       scaled_threshold = threshold * peak_gauss;
  107.       if (peak_gauss / sum < scaled_threshold)
  108.          scaled_threshold = peak_gauss / sum;
  109.       for (k=0; k<SQR(region); k++)           /* Nearby pixel centers */
  110.          {
  111.      if (kernel_tmp[k]/sum >= scaled_threshold)
  112.         inside_sum += kernel_tmp[k];
  113.      }
  114. /* Normalize coefficients exceeding threshold and fill in kernel arrays */     
  115.      
  116.       *ks++ = (long)(kc - kernel.coeffs);     /* x,y coeffs start this far down array */
  117.       k = 0;
  118.       for (y=0; y<region; y++)                /* Nearby pixel centers */
  119.       for (x=0; x<region; x++) 
  120.          {
  121.      if (kernel_tmp[k]/sum >= scaled_threshold)
  122.         {
  123.         *ko++ = x - region/2;             /* x,y offset from center pixel */
  124.         *ko++ = y - region/2;
  125.         *kc++ = kernel_tmp[k]/inside_sum; /* Normalized coefficient */
  126.         coeff_count++;
  127.         }
  128.      k++;
  129.      }
  130.       }
  131.    *ks = (long)(kc - kernel.coeffs);          /* One last one to mark end of final entry */
  132.    printf("%d x %d subpixel grid, %d x %d pixels considered.  Average of %5.1f coeff/sample\n",
  133.    sgrid,sgrid,region,region,(float)coeff_count/(float)SQR(sgrid));
  134.       
  135.    free (kernel_tmp);                         /* Done with this */
  136. }
  137.  
  138. /***************************************************************************/
  139.  
  140. long sample_gaussian(image_t *src, float px, float py)
  141. {
  142. long i,num_subpixels,list_id,list_offset,next_list_offset,count,src_pixel;
  143. long integer_x,subpixel_grid_offset_x,src_x,maxx;
  144. long integer_y,subpixel_grid_offset_y,src_y,maxy;
  145. float frac_x,frac_y,red_sum,grn_sum,blu_sum,red,grn,blu,coeff;
  146. float *coeff_ptr;
  147. short *offset_ptr;
  148.  
  149.    maxx = src->whole_xsize-1;
  150.    maxy = src->whole_ysize-1;
  151.    if (px < 0)    px = 0.0;
  152.    if (px > maxx) px = (float)maxx;
  153.    if (py < 0)    py = 0.0;
  154.    if (py > maxy) py = (float)maxy;
  155.    integer_x = (long)px;
  156.    integer_y = (long)py;
  157.    frac_x = px - (float)integer_x;    /* Between 0 and .9999999 */
  158.    frac_y = py - (float)integer_y;    /* Between 0 and .9999999 */
  159.    num_subpixels = kernel.subpixel_grid_size;
  160.    subpixel_grid_offset_x = (long)(frac_x * (float)num_subpixels);
  161.    subpixel_grid_offset_y = (long)(frac_y * (float)num_subpixels);
  162.    list_id = subpixel_grid_offset_x + subpixel_grid_offset_y * num_subpixels;
  163.    list_offset = *(kernel.starts + list_id);
  164.    next_list_offset = *(kernel.starts + list_id+1);
  165.    count = next_list_offset - list_offset;    /* This many coefficients */
  166.    coeff_ptr = kernel.coeffs + list_offset;
  167.    offset_ptr = kernel.offsets + 2 * list_offset;
  168.    
  169.    red_sum = 0.0;
  170.    grn_sum = 0.0;
  171.    blu_sum = 0.0;
  172.  
  173.    for (i=0; i<count; i++)
  174.       {
  175.       src_x = integer_x + (long)*offset_ptr++;
  176.       src_y = integer_y + (long)*offset_ptr++;
  177.       if (src_x < 0)    src_x = 0;
  178.       if (src_x > maxx) src_x = maxx;
  179.       if (src_y < 0)    src_y = 0;
  180.       if (src_y > maxy) src_y = maxy;
  181.       
  182.       coeff = *coeff_ptr++;
  183.       src_pixel = *IMG_INDEX(src,src_x,src_y);
  184.       red = (float)EXTRACT_RED(src_pixel);
  185.       grn = (float)EXTRACT_GRN(src_pixel);
  186.       blu = (float)EXTRACT_BLU(src_pixel);
  187.       
  188.       red_sum += red * coeff;
  189.       grn_sum += grn * coeff;
  190.       blu_sum += blu * coeff;
  191.       }
  192.       
  193.    return(MAKE_LARGE_RGB((long)red_sum,(long)grn_sum,(long)blu_sum));
  194. }
  195.  
  196. /***************************************************************************/
  197.  
  198. /* Copy image with arbitrary warp and/or zoom.  Use filter when sampling   */
  199. /* source image.  Traverse qp quad linearly to address source image.       */
  200. /* Traverse allocated memory to address dest image.  Dest image will be    */
  201. /* dest_width x dest_height pixels.  Outer 2-pixel border will be garbage, */
  202. /* "Undo_image" will be source, "Improved_image" destination.              */
  203.  
  204.       
  205. void *resample_quad(quad_points *qp, long dest_width, long dest_height)
  206. {
  207. long i,j,k,x,y,larger_dim,first_row,last_row,first_col,last_col,num_rows,num_cols;
  208. long *src_ptr,*dest_ptr,*dbg;
  209. long dest_border_x,dest_border_y;
  210. float fwb,fhb,xa,xb,xc,xd,xe,xf,xg,xh,ya,yb,yc,yd,ye,yf,yg,yh;
  211. float fw,xl,xr,xi,ldx,rdx,dx;
  212. float fh,yl,yr,yi,ldy,rdy,dy;
  213. float xzoom,yzoom,zoom,qp_rescale,minx,miny,maxx,maxy;
  214. int d;
  215.  
  216.    if (no_crop)
  217.       {
  218.       dest_border_x = 0;
  219.       dest_border_y = 0;
  220.       if ((dest_width != improved_image->whole_xsize) || (dest_height != improved_image->whole_ysize))
  221.          {
  222.          free_image(improved_image);
  223.          improved_image = make_new_image(dest_width,dest_height);
  224.          }
  225.       }
  226.    else
  227.       {
  228.       dest_border_x = dest_width / 14;
  229.       dest_border_y = dest_height / 14;
  230.       if ((dest_width != improved_image->interior_xsize) || (dest_height != improved_image->interior_ysize))
  231.          {
  232.          free_image(improved_image);
  233.          improved_image = make_new_image(dest_width + 2 * dest_border_x,dest_height + 2 * dest_border_y);
  234.          }
  235.       }
  236.    
  237. /* qp describes source quadrilateral in screen coordinates. The source    */
  238. /* image (undo_image) has arbitrary size, but it was resampled so largest */
  239. /* dimension was 512 when it was displayed (i.e., rescaled by "zoom").    */
  240. /* This rescaled image retains the original aspect ratio.  To get qp to   */
  241. /* align with the source image, we must subtract the window image origin  */
  242. /* and scale it so 0-511 on screen is 0-largest_dimension in source image.*/
  243.  
  244.    larger_dim = MAX2(undo_image->interior_xsize,undo_image->interior_ysize);
  245.    qp_rescale = (float)larger_dim / 512.;
  246.    minx = miny = 999;
  247.    maxx = maxy = -999;
  248.    for (i=0; i<8; i++)
  249.       {
  250.       (*qp)[i] -= IMAGE_BORDER;   
  251.       (*qp)[i] *= qp_rescale;
  252.       (*qp)[i] -= 0.5;   /* Point-sampled image has integer pixels */
  253.       if (EVEN(i))
  254.          {
  255.          (*qp)[i] += undo_image->x_border;
  256.          minx = MIN2(minx,(*qp)[i]);
  257.          maxx = MAX2(maxx,(*qp)[i]);
  258.      }
  259.       else
  260.          {
  261.          (*qp)[i] += undo_image->y_border;
  262.          miny = MIN2(miny,(*qp)[i]);
  263.          maxy = MAX2(maxy,(*qp)[i]);
  264.      }
  265.       }
  266.    
  267. /* If "crop_1_to_1" is asserted, then no resampling is to be done...just copy */   
  268.    
  269.    if (crop_1_to_1)
  270.       {
  271.       first_row = (long)((*qp)[1] - (float)dest_border_y + 0.5);
  272.       last_row =  (long)((*qp)[5] + (float)dest_border_y - 0.5);
  273.       first_col = (long)((*qp)[0] - (float)dest_border_x + 0.5);
  274.       last_col =  (long)((*qp)[4] + (float)dest_border_x - 0.5);
  275.       src_ptr = undo_image->whole_org + first_col + first_row * undo_image->pixels_per_line;
  276.       dest_ptr = improved_image->whole_org;
  277.       num_rows = ABS(last_row - first_row + 1);
  278.       num_cols = ABS(last_col - first_col + 1);
  279.       for (y=first_row; y<=last_row; y++)
  280.          {
  281.      memcpy(dest_ptr,src_ptr,num_cols * sizeof(long),dest_ptr - improved_image->whole_org);
  282.      src_ptr += undo_image->pixels_per_line;
  283.      dest_ptr += improved_image->pixels_per_line;
  284.      percentdone(y - first_row, num_rows);
  285.      }
  286.       return;
  287.       }
  288.       
  289. /*   qp subscripts: 6,7---4,5              */      
  290. /*                   / |   \               */
  291. /*                  / |     \              */ 
  292. /*               0,1---------2,3           */
  293. /*                                         */
  294. /*                                         */ 
  295. /*   Border labels:   d         h          */
  296. /*                                         */
  297. /*                  c   67---45   g        */
  298. /*                      /     \            */
  299. /*                     /       \           */
  300. /*               b   01---------23   f     */
  301. /*                                         */
  302. /*             a                       e   */
  303.  
  304. /*  qp is source image coordinates selecting quadrilateral which becomes destination interior */
  305. /*  a/e/h/d is source image coordinates which become destination whole corner pixels */
  306. /* They are proj_border away from qp, where proj_border is the projection of the destination */
  307. /* border into the source image */
  308.  
  309.    
  310.    xzoom = (float)dest_width / (maxx - minx);
  311.    yzoom = (float)dest_height / (maxy - miny);
  312.    zoom = AVG2F(xzoom,yzoom);
  313.    
  314.    generate_kernels(0.4/zoom,0.005);             /* Generate lists of coefficients */
  315.    
  316.    fwb = (float)(dest_width - 1 + 2 * dest_border_x);   /* Whole image width */
  317.    fhb = (float)(dest_height - 1 + 2 * dest_border_y);  /* Whole image height */
  318.    
  319. /* Extrapolate along 01 <--> 23 and 67 <--> 45 to get points b, c and f, g */
  320. /* Assume 1/16th border around rubber band "qp" is the source image border */
  321.  
  322.    if (no_crop)
  323.       {
  324.       xa = (*qp)[0];
  325.       ya = (*qp)[1];
  326.       xe = (*qp)[2];
  327.       ye = (*qp)[3];
  328.       xh = (*qp)[4];
  329.       yh = (*qp)[5];
  330.       xd = (*qp)[6];
  331.       yd = (*qp)[7];
  332.       }
  333.    else
  334.       {
  335.       xc = (*qp)[6] - ((*qp)[4] - (*qp)[6]) / 14.0;
  336.       yc = (*qp)[7] - ((*qp)[5] - (*qp)[7]) / 14.0;
  337.       xb = (*qp)[0] - ((*qp)[2] - (*qp)[0]) / 14.0;
  338.       yb = (*qp)[1] - ((*qp)[3] - (*qp)[1]) / 14.0;
  339.       xg = (*qp)[4] + ((*qp)[4] - (*qp)[6]) / 14.0;
  340.       yg = (*qp)[5] + ((*qp)[5] - (*qp)[7]) / 14.0;
  341.       xf = (*qp)[2] + ((*qp)[2] - (*qp)[0]) / 14.0;
  342.       yf = (*qp)[3] + ((*qp)[3] - (*qp)[1]) / 14.0;
  343.    
  344. /* Extrapolate along b <--> c and f <--> g to get points a, d and e, h */
  345.   
  346.       xd = xc + (xc - xb) / 14.0;   
  347.       yd = yc + (yc - yb) / 14.0;
  348.       xa = xb - (xc - xb) / 14.0;
  349.       ya = yb - (yc - yb) / 14.0;
  350.       xh = xg + (xg - xf) / 14.0;
  351.       yh = yg + (yg - yf) / 14.0;
  352.       xe = xf - (xg - xf) / 14.0;
  353.       ye = yf - (yg - yf) / 14.0;
  354.       }
  355.    
  356. /* We are going to traverse source region in a,e,d,h order */   
  357. /* "i" is source starting point, point a                   */
  358. /* Left side "l" starts there also                         */
  359. /* Right side "r" starts at point e                        */
  360.       
  361.    xi = xl = xa;
  362.    yi = yl = ya;
  363.    xr = xe;
  364.    yr = ye;
  365.    
  366. /* Rate of change along left edge is  (d-a)/fhb */   
  367. /* Rate of change along right edge is (h-e)/fhb */   
  368.    
  369.    ldx = (xd - xa) / fhb;
  370.    ldy = (yd - ya) / fhb;
  371.    rdx = (xh - xe) / fhb;
  372.    rdy = (yh - ye) / fhb;
  373.    
  374.    dest_ptr = improved_image->whole_org;     /* Storage for the cropped image */
  375.    
  376.    for (y=0; y<dest_height + 2 * dest_border_y; y++)
  377.       {
  378.       dx = (xr - xl) / fwb;                  /* Compute new scanline increments */
  379.       dy = (yr - yl) / fwb;
  380.       for (x=0; x<dest_width + 2 * dest_border_x; x++)
  381.          {                                   /* Sample & copy */
  382.          *dest_ptr++ = sample_gaussian(undo_image,xi,yi);
  383.          xi += dx;                           /* Bump scanline address */
  384.          yi += dy;
  385.      }
  386.       xi = xl += ldx;                        /* Advance edge walk addresses */
  387.       yi = yl += ldy;
  388.       xr += rdx;
  389.       yr += rdy;
  390.       percentdone(y, dest_height + 2 * dest_border_y);
  391.       }
  392.  
  393.    free(kernel.starts);
  394.    free(kernel.offsets);
  395.    free(kernel.coeffs);
  396. }      
  397.  
  398.  
  399. /********************************************************************************/
  400.